Dynamics in the Sherrington-Kirkpatrick Ising spin glass at and above T g 
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A detailed numerical study is made of relaxation at equilibrium in the Sherrington-Kirkpatrick 
Ising spin glass model, at and above the critical temperature T g . The data show a long time stretched 
exponential relaxation q(t) ~ exp [-(t/r(T))' 3(T) ] with an exponent j3(T) tending to « 1/3 at T g . 
The results are compared to those which were observed by Ogielski in the 3d ISG model, and are 
discussed in terms of a phase space percolation transition scenario. 
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INTRODUCTION 



The Sherrington-Kirkpatrick (SK) Ising spin glass 
(ISG) model has been intensively studied for almost forty 
years. Introduced [H, 0] as a starting point for study- 
ing Edwards-Anderson-like Q spin glasses, it is a clas- 
sical mean-field TV-spin Ising model in which there are 
quenched random interactions between all pairs of spins. 
It has an ordered spin glass phase below the critical tem- 
perature T g = 1. The static properties such as the spe- 
cific heat, the magnetic susceptibility and the spin glass 
order parameter are known to high precision in the or- 
dered state and are well understood in terms of the Parisi 
Replica Symmetry Breaking (RSB) theory [3] and its 
subsequent developments. The paramagnetic regime has 
been considered "trivial" in that there are simple exact 
expressions for the basic static physical properties in the 
large N limit. 

In the ordered state the dynamics are very slow; the 
equilibrium relaxation time diverges exponentially when 
the number of sites N goes to infinity - For the 
paramagnetic state, the dynamics were discussed in early 
work, Refs. @, [lMl- The ISG relaxation corre- 
sponds to a continuous spectrum of relaxation frequen- 
cies, for which an explicit linearized expression was given 
in Ref. Q, leading to a relaxation function for Gaussian 
interactions and Glauber dynamics in the form of an inte- 
gral, Ref. 0, Eq. 5.27. It was concluded that in the large 
size limit at T g the relaxation at thermal equilibrium of 
the autocorrelation q(t) = (1/-/V) Z)j < Si(0).Si(t) > as 
a function of time t (in Monte Carlo steps per spin) 
would be q(t) <~ i -1 / 2 , and above T g the very long 
time limit would take the form of a simple exponential 
q(t) ~ exp[— (t/r(T))] corresponding to a cutoff in the 
relaxation frequency spectrum. As far as we are aware, 
after exploratory numerical calculations on a very lim- 
ited number of Gaussian interaction samples over a few 
Glauber time steps in Ref. Q no further numerical work 
on the dynamics in this temperature region has been re- 
ported. 

Here a detailed numerical study of the thermal equi- 
librium relaxation of the autocorrelation function q(t, N) 



in the ordered state as a function of time t in updates 
per spin and N the number of spins is extended to 
temperatures at and above T g . At T g the relaxation is 
always hnite size limited; the expected scaling form 



g(i,A0A^ 1/3 = F[-(t/N 2 / 3 )] 



(1) 



is observed. For the three temperatures above T g at 
which measurements were made the large size limit be- 
havior [q(t, N) independent of N] is reached with the 
largest samples studied. 

For q(t) values down to at least q(t) ~ 10~ 6 the au- 
tocorrelation function decay is strongly non-exponential. 
Satisfactory long t fits are given both above T g and in the 
finite size scaling regime at T g by stretched exponentials 



q (t) ~ ex P [-(t/r(T))W] 



(2) 



with a temperature dependent exponent /3(T) which is 
smaller than 1 and which tends to j3(T) <~ 1/3 at T = T g 

M. 



THE MEAN FIELD FERROMAGNET 

Before discussing the SK data analysis it is instructive 
as an illustration to recall the (non-equilibrium) relax- 
ation behavior for the mean field Ising fcrromagnet, for 
which exact results by Suzuki and Kubo [l4[ were dis- 
cussed in Ref. 0]. At temperatures above T c = 1 the 
relaxation takes the form 

m{t,T) = 

(1 - /3) 1/2 /[(l -P + P 3 /3) exp(2(l - fit) - /? 3 /3] 1/2 



where /3 = l/T, which becomes 

m(t,T c ) = (l + (2/3)i)- 1 / 2 



(3) 



(4) 



at T c . It is easy to see from Eq. [3] that above T c and 
beyond a time t ~ t(T) = 1/(1-/3), m(t, T) is essentially 
equal to A(T) exp(— t/r(T)), a pure exponential with a 
time independent prefactor, 



A(T) = [(1 -/?)/(! -£ + /? 3 /3)] 1/2 



(5) 
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FIG. 1: (Color online) An example of an m(t) plot for the 
mean field ferromagnet calculated using Eq. [3] The points 
correspond to m(t) at T = f.05T c ; the black line is the alge- 
braic critical behavior, and the blue curve is the asymptotic 
pure exponential. 

see Figure 1. Beyond t ~ r only one single mode con- 
tributes to the relaxation in the mean field ferromagnet. 
At short times (time zero to time ~ t(T)) there is a more 
complex function which connects up between the starting 
point m(0, T) = 1 and the asymptotic region. At no time 
for temperatures above T c is m(t,T) « t~ x l 2 exp(— t/r) 
a useful approximation. 

It can be helpful to show the data in the form of 
a derivative plot — dlog(m(t)) / dlog(t) against t. Then 
a pure exponential m(t) = Acxp[— t/r] appears as a 
straight line of slope 1/r passing through the origin. 

SK SIMULATIONS 

We now return to the SK simulations. Systems of size 
N = 64, 128, 256, 512, 1024 and 2048 were equilibrated 
using the procedure described in Ref. Simulations 
were made on 1024 independent disorder samples at each 
size with two clones simulated in parallel. (As shown in 
Ref. Q two clones allow one to build up a noisy but 
unbiased estimate of the thermal fluctuations. Here they 
are only used to decrease the thermal fluctuations). We 
performed 10 7 Monte Carlo updates per site after ther- 
malization, with measurements of q(t) made every 4 time 
intervals. Both Metropolis and heat bath relaxation runs 
were carried out. It can be noted that SK simulations are 
numerically demanding : there are about as many indi- 
vidual spin-spin interactions in an N = 2048 SK sample 
as in an L = 90 sample in dimension 3. 

In the TV = 256 data set, and that set only, it turns out 



that one specific disorder sample i gives qi (t, T) that are 
very far from the values for all the 1023 other samples 
of the same size. This is true for sample i with either 
Heat Bath or Metropolis updating schemes, and for all 
values of T. For example at T = T g = 1 we have 1024 
estimates for g(256,£ = 1024) with a median equal to 
0.0000040, first and third quartiles equal to -0.0001183 
and 0.0001292 respectively, and an isolated maximum 
value equal to 0.2946040. This single outlier was omitted 
from the global analysis, since it would have led to a long 
t distortion of our values of the mean < <?(256, t, T) > and 
the statistical error. At long t when < q(t,T) > is very 
small, if one single sample by statistical accident is very 
uncharacteristic it can have a big influence on the mean 
and the error. Note that in our simulation the thermal 
noise on < q(t, T) > is much smaller than the disorder 
noise. 

The temperatures at which measurements were made 
are T = T g = 1 and T = 1.1, 1.2, 1.3 

SK ANALYSIS 

The standard dynamic finite size scaling rule exactly 
at the critical point can be written 

q c (L,t)=L-^F[t/r(L)} (6) 

with t(L) ~ L z . It is well known that quite generally 
finite size scaling is modified above the upper critical di- 
mension d u , with the combination (T — T g )L 1 l v replaced 
by (T - T^A^/K'') [§], where N is the number of sites, 
namely L is replaced by N x / du . 

The SK model is the infinite dimensional version of a 
model with upper critical dimension d u = 6 9J, and its 
scaling behavior is accordingly described by the effective 
length scale L e // = iV 1 / 6 . This result can be obtained 
directly by noting that in the critical region the Landau 
expansion of the free energy of the replicated model starts 
(symbolically) by a(T ~ T g )Q 2 + bQ 3 with Q the n by 
n matrix q a ^ where n is the number of replica. The 
partition function takes the form J dQexp(—N(a(T — 
T g )Q 2 + bQ 3 )), and after the change of variable Q = 
X/[(T — Tg)^] 1 / 2 one obtains a free energy that is a 
function of (T — T g )N 1 / 3 , in agreement with the above 
statement since v = 1/2 for this model. 

Therefore for dimension 6 and above (including SK) 
the effective ISG "length" scale is L eff = N 1 / 6 . The 
mean field ISG critical exponents are j3 = 1, v = 1/2, 
and z = 4. 

Hence in the SK case the dynamic scaling becomes 

q c (N,t) - N^^Fit/N 2 ^ 3 } (7) 

so that in a critical temperature scaling plot one should 
expect q c (N,t)Ny 3 to be a function F[t/N 2 / 3 ]. 
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FIG. 2: (Color online) The short time scaled SK relax- 
ation data q^N 1 / 3 against t/N 2 ^ 3 in logarithmic coordi- 
nates at the critical temperature T = 1. The colors for 
N = 64, 128, 256, 512, 1024, 2048 are yellow, pink, cyan, green, 
red and black. For clarity the iV = 2048 points are larger than 
the others. The blue line indicates q c (t) ~ £ -1//2 . 



FIG. 3: (Color online) The scaled SK relaxation data q(t)N 1/3 
in logarithmic coordinates against t/N 2 ^ 3 at the critical tem- 
perature T = 1. The color code is the same as in Figure 
2. The blue curve is a stretched exponential q{t)N 1 ^ 3 = 
14exp[-((t/iV 2/3 ) 1/3 )]. 



In conditions where q c (N,t) approaches the N — > oo 
limit, the scaling rule becomes 

q c (t) ~ N-^it/N 2 ' 3 )- 1 ' 2 = r 1 ' 2 (8) 

independent of N. This is the situation for short times 
and large N, Figure 2. 

For reasons to be discussed below, in addition to test- 
ing this scaling form we want to test if the long time 
critical finite size scaling function F[x) is consistent with 
the stretched exponential form B exp[— (x)' 3 ] having the 
particular exponent value ft = 1/3. In figure 3 the 
scaled critical temperature relaxation data arc plotted 
as log^^iV 1 / 3 ) on the y-axis against (t/N 2 ^ 3 ) on the 
x-axis. The long time critical scaling rule is well obeyed 
over the whole range of sizes N which have been studied; 
beyond times of a few Monte Carlo steps q c {N, t) TV 1 / 3 is 
indeed an iV-indcpendent function F[t/N 2 / 3 ] to within 
the error bars for all t and N, with no visible sign of 
corrections to finite size scaling. 

The initial very short time relaxation (Figure 2, invisi- 
ble in Figure 3) follows approximately the iV-independent 
form q(t) ~ t -1 / 2 as predicted. As can be seen in Figure 
3, within the numerical precision the finite size scaling 
function is from then on compatible with a stretched ex- 
ponential having exponent (3 = 1/3. 

For the three temperatures above T g , T = 1.1, 1.2, 1.3 
the raw q(N, t) data are shown in Figures 4, 5 and 6. 
It can be seen that at each of these temperatures the 
size dependence saturates with increasing N. The q(N, t) 
curves for the largest sizes coincide within the numerical 



error bars, so these data can be taken to represent the 
"infinite" N limit behavior. 

The q(t) curves calculated with the expression given 
in Ref. [2| Eq. 5.27, are almost identical to the limiting 
numerical q(N,t) curves for T = 1.1 and T = 1.2 but 
the calculated curve for T = 1.3 lies somewhat above the 
numerical data and has a rather different shape. The 
good quantitative agreement for the first two tempera- 
tures is perhaps fortuitous as the expression in Ref. Q 
corresponds to a model with Gaussian interactions, and 
a linearization approximation was made in the calcula- 
tion. The authors expected the expression to become 
accurate only at temperatures well above T g . Neverthe- 
less the approach of Ref. is substantially validated by 
the numerical data. 

We can note that the simulations made with the 
Metropolis updating rules (not shown) provide q(t, N) 
data which are very similar to the heatbath data but 
with a global shift in the time scale to times shorter by a 
factor of between 2 and 3 depending on the temperature. 
It can be shown that as a general rule q(t) will fall faster 
with Metropolis updating than with heatbath/ Glauber 
updating; the Metropolis acceptance rate is higher than 
that for heatbath/Glauber with a ratio of which depends 
on temperature. 

The figures are of the same form as in the ferromag- 
netic example of Figure 1 with an initial short time 
regime followed by a long time asymptote, except that 
the long time fcrromagnet pure exponential asymptote is 
replaced by a stretched exponential. Beyond the initial 
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FIG. 4: (Color online) The raw relaxation data q(t) against t 
at the temperature T — 1.1. The color code for N is the same 
as in Figure 2. The fit Eq. [2] to the largest sizes is shown as 
the full black curve. q(t) calculated from Ref. 0| Eq. 5.27 is 
shown as the full blue curve. 
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FIG. 5: (Color online) The raw relaxation data q(t) against t 
at the temperature T — 1.2. The color code for N is the same 
as in Figure 2. The fit Eq. [2] to the largest sizes is shown as 
the full black curve. q(t) calculated from Ref. 0| Eq. 5.27 is 
shown as the full blue curve. 

short time regime the stretched exponentials fit the large 
N SK data to within the numerical precision. 

If the stretched exponential form Eq. [2] is assumed, 
then 

- dlog(?(t))/dbg(t) = P{T)(t/rf {T) (9) 
An example of a plot of — d\og{q(t)) / d\og(t) against t 
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FIG. 6: (Color online) The raw relaxation data q(t) against t 
at the temperature T — 1.3. The color code for N is the same 
as in Figure 2. The fit Eq. [2] to the largest sizes is shown as 
the full black curve. q(t) calculated from Ref. Q| Eq. 5.27 is 
shown as the full blue curve. 
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FIG. 7: (Color online) The derivative plot 

— dlog(q(t)) / d\og(t) against t at T = 1.1 for the two 
largest sizes, N = 2048 (black) and 1024 (red). The fit (blue 
curve) is for the same parameters as in Figure 4. 

is shown in Figure 7. For a pure exponential relax- 
ation j3 = 1 the data points would lie on a straight 
line going through the origin. This is clearly not the 
case. By analogy with the ferromagnetic behavior, be- 
yond an initial short time region the large N data are 
fitted by a B(T)t^ curves with (3{T) less than 1 and 
B(T) ^ I3{T)/t{T)^ t \ 
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Satisfactory fits are obtained with : /3(1.1) « 0.52, 
r(l.l) w 6.8, 13(1.2) « 0.61, r(1.2) « 4.2,and £(1.3) w 
0.68, r(1.3) sw 3.2 for the heatbath updates. These fits 
are shown as the black curves in Figures 4 to 6. 

The relaxation time t(T) is of course increasing as T g is 
approached but with only three temperature points it is 
difficult to estimate a functional form for t(T). However 
the most significant result is that the relaxation above T g 
in the effectively infinite size limit can be represented by 
stretched exponentials having exponents (3(T) which are 
less than 1 and which drop regularly as T g is approached. 

From this result for the large N limit relaxation data 
at the temperatures above T g together with the scaled 
finite size limited relaxation at criticality it can thus 
be concluded that in the paramagnetic region the SK 
model thermodynamic limit equilibrium relaxation can 
be treated as essentially stretched exponential, with an 
exponent (3(T) which decreases continuously as the tem- 
perature is lowered. (3(T) tends to a value /3 C w 1/3 
at criticality. It should be underlined that this behavior 
has been followed down numerically though five orders of 
magnitude in q(t) at each temperature. 



ISG IN DIMENSION THREE 

It is instructive to compare with the situation for a 
finite dimension model. Ogielski [l5| fitted thermody- 
namic limit q(t) data on the bimodal ISG in dimension 
3 using an empirical function of the form 



q(t) ~ t-*W exp[-(i/r(T))' 3 ( T )] 



(10) 



where an algebraic prcfactor multiplies a stretched expo- 
nential. Ogielski estimated Tg « 1.175 which is close to 



recent estimates T g m 1.12 3li7j- The values estimated 



for the stretched exponential exponent /3(T) decrease reg- 
ularly from f3(T) « 1 at a temperature of about 4.0 to 
/3(T) ~ 0.35 at a temperature just above T g . The relax- 
ation time t(T) diverges as T — > T g and the algebraic 
prefactor exponent x(T) decreases from « 0.5 at high 
temperatures to (d - 2 + rj)/2z ~ 0.070 at T g . These 
equilibrium relaxation results were obtained from high 
quality simulations on very large samples (up to L = 64) 
and have never been improved on since. 

Thus the general pattern of the paramagnetic state re- 
laxation reported by Ogielski in the 3d ISG and the ob- 
served behavior of the SK model described above show 
striking similarities; above all both sets of data in the 
paramagnetic region are consistent with stretched expo- 
nential decay having an exponent (3(T) tending to a crit- 
ical value (3 C ~ 1/3 at T g . These two models represent 
the lowest integer dimension at which finite temperature 
ISG ordering takes place and the infinite dimension mean 
field limit (for which the RSB theory is well established) 
respectively, so it would seem reasonable to expect that 



the same pattern of relaxation should hold for ISGs at 
all intermediate dimensions also, both above and below 
the upper critical dimension (see Ref. [HI). 

Shortly after Ogielski's work was published it was 
noted that his data were consistent with (3(T g ) = 1/3, 
and it was conjectured (lij that stretched exponential 
relaxation with an exponent tending to precisely 1/3 at 
criticality could be universal in ISGs. The argument is 
briefly summarized in the next section. 



PERCOLATION TRANSITION SCENARIO 

A thermodynamic phase transition can be regarded as 
a qualitative change in the topology of the thermody- 
namically attainable phase space with decreasing tem- 
perature, i.e. with decreasing internal energy. Thus for 
a standard ferromagnetic transition a high temperature 
spherical phase space becomes more and more "elliptic" 
as the temperature is lowered and the number of attain- 
able states drops; finally at T c it splits into two (up and 
down) mirror-image subspaces. In a naive scenario based 
on RSB, for an ISG at T g phase space shatters into a large 
number of inequivalent clusters. In Euclidean space there 
is a well studied transition of this type, the percolation 
transition. For a concentration p > p c there is a giant 
cluster of sites while just below p c there are only small 
inequivalent clusters . 

The total phase space of an TV-spin S = 1/2 Ising 
system is an TV dimensional hypercube. Relaxation of 
any TV-spin Ising system by successive single spin up- 
dates can be considered strictly as a random walk of 
the system point along near neighbor edges among the 
thermodynamically attainable vertices on this hypercube 
[TBI . It was argued that as random walks on full 
[hyper] s phe rical surfaces result in pure exponential de- 



cay [20|, |2l| and random walks on threshold percolation 



clusters in Euclidean space lead to sub-linear diffusion 
(R 2 ) ~ t^ d [22|], random walks on threshold percola- 
tion clusters inscribed on [hyper] spheres would be char- 
acterized by "sub-exponential" relaxation of the form 
q(t) = (cos(0(t))) = exp[— {t/r) l3d ] with the same expo- 
nents /3d as in the corresponding Euclidean space. This 
was demonstrated numerically for d = 3 to 8 [23j . A hy- 
percube being topologically equivalent to a hyperspherc, 



on a diluted hypercube at threshold [24j, |25j random 
walks can be expected to lead to stretched exponential 
relaxation with exponent /3 = 1/3. Successive explicit 
numerical studies of random walks on the randomly oc- 
cupied hypercube at p c [26T - [28| have confirmed this. 

It is important to note that the limiting behavior for re- 
laxation due to diffusion on sparse graphs has been shown 
analytically to take the form of a stretched exponential 
with exponent /3 = 1/3 0, l3p|. 

It was further conjectured [191 ] that in a complex Ising 
system such as an ISG the phase transition would be 
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analogous to a percolation transition in configuration 
space. Detailed "rough landscape" models for the con- 
figuration space of complex systems have been widely in- 
voked (see for instance Refs. [3l|-[33j|); these models can 
be thought of in terms of linked basins with a gradual di- 
lution of the links leading finally to a percolation thresh- 
old. As the critical properties of a percolation transi- 
tion are very robust, if the configuration space perco- 
lation threshold scenario is valid the basic critical be- 
havior should not be sensitive to model details; in par- 
ticular stretched exponential relaxation with exponent 
1/3 should be observed generically in ISGs whatever the 
space dimension and possibly also in a wider class of com- 
plex systems. 



CONCLUSION 

The present numerical SK results taken together with 
Ogielski's 3d ISG analysis [HI can be taken as a strong 
empirical indication of a universal equilibrium relaxation 
pattern for ISGs in the paramagnetic regime : stretched 
exponential decay Eq. [2] having an exponent j3(T) which 
tends to 1 /3 when the relaxation time diverges at the crit- 
ical temperature T g . Auto-correlation function decay at 
the percolation transition on the randomly occupied hy- 
percube [26-28 has precisely this form, suggesting that 
the ISG phase transition can indeed be considered in 
terms of a percolation transition in phase space. 

The same type of relaxation pattern has been observed 
in many other complex system studies, both numerical 
and experimental. It is tempting to conclude that the 
physical scenario could have implications for a wide class 
of transitions, extending well beyond the ISG family. 
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